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ABSTRACT 

We investigate the long-time evolution of magnetic fields generated by the two-stream instability at 
ultra- and sub-relativistic astrophysical collisionless shocks. Based on 3D PIC simulation results, we 
introduce a 2D toy model of interacting current filaments. Within the framework of this model, we 
demonstrate that the field correlation scale in the region far downstream the shock grows nearly as 
the light crossing time, A(t) ~ ct, thus making the diffusive field dissipation inefficient. The obtained 
theoretical scaling is tested using numerical PIC simulations. This result extends our understanding of 
the structure of collisionless shocks in gamma-ray bursts and other astrophysical objects. 

Subject headings: shock waves — magnetic fields — gamma rays: bursts — supernova remnants 
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1. INTRODUCTION 

Internal and external shocks in gamma-ray bursters 
(GRBs), internal shocks produced in jets of micro- and 
normal quasars and in active galactic nuclei jets, shocks 
in supernovae remnants, merger shocks in galaxy clusters 
and large scale structure, — all of them represent a single 
class of strong collisionless shocks. The theoretical predic- 
tion that the counterstreaming instabili ty operating at the 
shock produces strong magnetic fields ijMedvedev &: Loea 
1999) has recently been confirmed in a number 
of st a te-of-the-art numerical s umulations (jSi lva. et all 
2003 1 iFrederiksen. et"afl l2004 iNishikawa. et all 12003: 
Saito fc Sakail 12004 iKazimura. et alJll998|) . In this pa- 
per we investigate the long-time, nonlinear evolution of 
the produced fields. The knowledge of the field dynam- 
ics downstream the shock is crucial for our understanding 
of the shock particle acceleration, as well as the spectral 
properties (e.g., variability) of radiation produced at as- 
trophysical shocks. 

2. THE MODEL 

Fully 3D PIC numerical simulations of the shocks 
demonstrate that the generated magnetic fields are associ- 
ated with a quasi-two-dime nsional distribution of current 
filaments ijSilva. et alJl2003|) . Hence we suggest the follow- 
ing toy model. 

We consider straight one-dimensional current filaments 
oriented in the vertical, z, direction. Initially, all filaments 
are identical: the initial diameter of them is Do, their ini- 
tial mass per unit length is fio ~ mn(itD\/ A), where m is 
the mass of plasma particles (e.g., electrons) and n is their 
number density. Each filament carries current Jo in either 
positive or negative z-direction. The net current in the 
system is set to zero, i.e., there are equal numbers of pos- 
itive and negative current filaments. We also assume that 
the initial separation (the distance between the centers, 
see Fig. ^) of the filaments do is comparable to their size, 
do ~ 2Dq. Finally, no external homogeneous magnetic 



field is present in the system. 

2.1. Single filament dynamics 

Initially, the filaments are at rest and their positions in 
space are random. This configuration is unstable because 
opposite currents repel each other, whereas like currents 
are attracted to each other and tend to coalesce and form 
larger current filaments. The characteristic scale of the 
magnetic field will accordingly increase with time. We 
study this process quantitatively using the toy model of 
two interacting filaments. There are two limiting cases: 
(i) when the characteristic velocity of the filaments during 
the coalescence process is much smaller than the speed of 
light and (ii) when these velocities are comparable. We 
consider these cases separately. 

2.1.1. Nonrelativistic motion 

The magnetic field produced by a straight filament is 
Bo(r) = 2Jo/(cr), where r is the cylindrical radius. The 
force per unit length acting on the second filament is 
dF/dl = —BqIq/c Since dF/dl — fix, where x is the 
position in the center of mass frame and "overdot" de- 
notes time derivative, we write the equation of motion as 
follows: 

oj 2 1 

* = (1) 

where we used that r = 2x and the reduced mass fi r = 
Ho/ 2. We define the coalescence time as the time re- 
quired for the filaments, which are initially at rest, to cross 
the distance between them and "touch" each other, which 
happens when the distance between their centers becomes 
equal to Do, i-e., when x — Dq/2. The coalescence time, 
as it is defined above, is independent of the details of the 
merging process itself, which involves rather complicated 
dynamics associated with the redistribution of currents. 
Quite obviously, the interaction between the filaments is 
the weakest at large distances x ~ xo ~ do/2. Hence, the 
coalescence rate is limited by the filament motions at the 
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largest scales. The coalescence time can be readily esti- 
mated from Eq. l|TJl. assuming that x ~ Xq ~ do/2 and 
i ~ (e?o/2)TQ~ 2 , as follows: 



t ,nr ~ (-DgcVo/(2/ 2 )) 



1/2 



(2) 



The above estimate is valid as long as the motion is non- 
relativistic. The maximum velocity of a filament is at the 
time of coalescence, x = Dq/2: 



v m o ~ Aj/2t ~ io/(c-v/2/zo)- 



(3) 



It must always be much smaller than the speed of light. 

It turns out that Eq. QJ can be solved exactly in terms 
of time as a function of the filament separation: 

t(x) = ^/jrjj^cxo/ (2J ) erte^y/-\n(x/x )) , (4) 

where xo = x(0) and erfc(y) = 2ir~ 1 / 2 dr\ is the 

complimentary error function. The coalescence time can 
be calculated exactly: 



to,nr = V^o c - D o/(4/o)erfc^Vln2j « 0.11^/^cD //o, 

(5) 

which is shorter than our estimate, Eq. (0, by a factor of 
seven. The velocity as a function of position is 



v(x) = 2/ /(cVmo) \/~ln(x/x ). 
Consequently, the maximum velocity is 



(6) 



v m0 = 2I /(cy 



1.67J /(<V£o)- (7) 



2.1.2. Relativistic motion 

If the motion of a filament during the merger becomes 
relativistic, the separation cannot decrease faster than as 
t(x) ~ x/c. Therefore, the coalescence time will be 



tq.r ^ (d /2)/c = D /c. 
2.2. Collective dynamics of filaments 



(8) 



We now consider the filament coalescence as a hierarchi- 
cal process. Suppose that initially the system contains Nq 
current filaments, with an average separation do ~ 2Z?o- 
Each of the filaments carries current Iq, its diameter is 
Do and its mass per unit length is fj,o- For simplicity, we 
asssume that filaments coalesce pair wise. 

Having the original "zeroth generation" of filaments 
merged (the process takes about tq,nr or tq,r to com- 
plete), the system will now contain No/2 of "first gener- 
ation" filaments. Each of these filaments carries current 
I\ = 2Io, has mass per unit length [i\ = 2fj,o, and the sepa- 
ration between them is d\ = v2^o (because the two dimen- 
sional number density of filaments decreased by 2). Since 
\i cx D 2 , the filament size also increases as D\ = \/2Do- 
Remarkably, this new configuration is identical to the ini- 
tial one, but with the re-scaled parameters. Hence, the 
coalescence process is self-similar. The produced first gen- 
eration filaments will be interacting with each other and 
merge again to yield the second generation. The coales- 
cence process will then continue in a self-similar way. Note 



that the coalescence times at each stage are not necessarily 
the same. 

At the k-th merger level, i.e., after k pairwise mergers, 
the filament current, its mass per unit length and its size 
are 



I k = 2 k Io, Mfc = 2Vo, D k = 2 fe / 2 A) 



(9) 



whereas the filament separation is dk ~ -Dfc/2. The coa- 
lescence time at fc-th level may be estimated in the way 
described in the prevous section. Using Eq. J2J or J3J and 
Eq. |JSJ|, we obtain 



Tk,NR — To t NRt T kM — 2 k l 2 To.R- 



(10) 



Since the coalescence time is independent of k while the fil- 
aments are non-relativistic, whereas the distance between 
them increases, the typical velocities of the merging fil- 
aments grow with time and, will approach c. From Eq. 
iJTJl, we obtain v m ^k = 2 k l 2 v m o- The transition from the 
non-relativistic regime to the relativistic one occurs after 
about k* mergers: 



21og 2 (c/u m0 ), 



(11) 



where v m o is set by the initial state of the system, Eq. Q . 

Our primary interest is the evolution of the transverse 
correlation length of the magnetic field. With the fil- 
amcts being randomly distributed in space, the character- 
istic scale on which the magnetic field fluctuates is about 
half the distance between the filaments, i.e., \s,k — dk/2. 
Thus, 



A s , fc ~ 2 k ' 2 Do. 



(12) 



Finally, it is instructive to present the evolution of the 
parameters as a function of physical time, t, rather than 
the merger level, k. Apparently, it takes t = J2k>=o T k' to 
complete k mergers, where Tk is given by Eq. IjlUI) . This 
equation implicitly defines k(t). In the non-relativistic 
case, the solution to this equation is obvious: k — t/ro, 
because tq is independent of k. In the relativistic case, the 
sum of the series is easily calculated using that Y^=o a * = 
(a n+1 — l)/(a — 1). Thus, for the non-relativistic and rel- 
ativistic cases respectively, the solution to the equation 
t = J2k'=o T k' reads 



k = t/To,NR, 

fc = 21og 2 \(t/T 0<R )(V2-l) + l 



(13) 
(14) 



Thus, the magnetic field correlation length increases as a 
function of time as 

<Vb(£) — Do2 t ^ 2r °- NR \ X B (t)^ct, (15) 

in the non-relativistic and relativistic regimes, respec- 
tively. Note that the last expression is an approximation 
at large times t ^> t ,r, i.e., at large k 3> 1. 
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3. APPLICATION TO COLLISIONLESS SHOCKS 

The growth rate and the saturation level of the field 
generated at shocks depend on the composition of the out- 
flowing: ionized gas llMedvedev fc Loeblll999t ISilva. et"all 
2003l lFrederiksen. et al.l2004jl . The gas composition in the 
cosmological outflows is not generally known. Considering 
the interaction of the ejecta with the interstellar medium, 
the external shock, it is quite reasonable to assume 
that the shock is propagating through an electron-proton 
plasma. In contrast, the ejecta itself, where internal shocks 
occur, may be either dominated by electron-positron pairs 
(leptonic jet) or by electrons and protons (baryonic jet) 
with or without e~e + -pair sub-population. 

Here, we will consider a general case of the ejecta con- 
taining several species, labeled by the subscript s. We as- 
sume that the species have different masses m s , but their 
charges, by the absolute value, are equal to e. In general, 
each species has the bulk Lorentz factor (in the center of 
mass frame) 7 S and the thermal Lorentz factor 7 S , the 
latter represents the random velocites of the particles (ini- 
tially, % < 7 S ). When the instability shuts off, the parti- 
cles are randomized over the pitch angle, hence 7 S ~ 7 S . In 
the context of astrophysical outflows, we assume that all 
the species have the same 7 S = T and 7 S ~ T, where T is 
either the Lorentz factor of the external shock in observer's 
frame, or the relative Lorentz factor of the two colliding 
shells measured in the frame comoving with the ejecta. In 
the latter case, the ejecta itself moves relativistically in 
observer's frame with the Lorentz factor r c j. 

It is convenient to introduce the equipartition parameter 
for each species: 

e B , s = B 2 /(87rn s j s m s c 2 ), (16) 

where B s is the strength of the magnetic field produced by 
the instability operating on the species s, n s is the num- 
ber density of particles (both parameters are measured in 
the shock comoving frame). The equipartition parameter 
describes the efficiency of the magnetic field generation 
process, i.e., what fraction of the total kinetic energy of 
the particles of each species goes into the magnetic field 
energy. Note that this definition of the equipartition pa- 
rameter is different from the commonly used definition, 
which describes the fraction of the total kinetic energy of 
the ejecta (summed over all species) that goes into mag- 
netic field. We define the plasma and Larmor frequencies 
as 

1/2 

ujp^s — (4:ire 2 n s /j s m s ) , (17) 
lu b ,s — eB s /j s m s c. (18) 

For convenience, we also introduce their "nonrelativistic" 

counterparts ui p>s = uj PtS ^ 2 and Cjb,s = ^B.sls- It is 
useful to remember the following relation: 

U B ,s/Up,s = \f^B~s- (19) 

We may now represent the main results of Section in 
terms of plasma parameters. First, the initial separation 
between the filaments, Dq, must be comparable to the 
characteristic correlation length of the magnetic field pro- 
duced by the instability. This length at the onset of the 
instability is, in turn, set by the wavenumber of the fastest 



growing mode ijMedvedev fc Loeblll999|) : A s — "fl^ 2 c/uj p _ s , 
which is a factor of (%/j s ) < 1 smaller than the relativis- 
tic skin depth, c/lu p , s . However, we cannot set this scale A s 
as Dp because, as indicated by 3D PIC e~e + simulations 
llSilva. et all2 003h €b is not constant in time at the begin- 
ning of the nonlinear filament interaction (at t ~ lOw"*)- 
Hence the analysis of Section [3 is not applicable in such a 
regime. In fact, it takes few plasma times, ui~],, in their 
simulations for cb to attain its asymptotic value. The field 
correlation scale at this moment (t ~ 15uj~1) is somewhat 
larger than A s . We include this uncertainty via the pa- 
rameter 77 > 1 as follows: 

D Q = rj(c/u> PtS ) = r)(cVf/u PjS ). (20) 

Second, using § B • d\ = (Att/c)I, we express the current 
Iq in terms of the equipartition parameter, eg s , as 

Jo =ryr^/e B)S /2(m s c 3 /e). (21) 

Third, the mass per unit length of a filament must take 
into account that the particle thermal motion is relativis- 
tic, hence the masses are 7 S ™ S ~ rm s . We obtain 

Ha = Tm s n s (ttD^/4) . (22) 

The temporal evolution of the field correlation scale is 
determined by Eq. where tq^r is given by iJSJ. The 

coalescence time may be written as 

a T t] _! a T r] _ a _ x 

to,nr = —z- u B = uj s ~uj s , (23) 

^ ^/XeB,s 

where we i ntrod uced the numerical factor a T = 
(•y/7r/4) erfc(Vha 2) » 0.11. Hereafter we assumed the typ- 
ical values: ts,s ~ 10~ 3 and 77 ~ l.The maximum merger 
velocity v m o, Eq. J7J), in terms of the speed of light is 

Vmo/c = a vy /8e B ,s ~ 0.1 , (24) 

wh ere another numerical factor is introduced: a v = 
2\/ln 2 w 1.67. The transition from the non-relativistic 
to relativistic coalescence regime occurs after fc* mergers, 
given by Eq. Ijlll) . i.e., at the time 

i* = 21og 2 (c/u 0m ) t , nr - 5t 0<N r. (25) 
4. NUMERICAL RESULTS 

We now compare our theoretical predictions with the 
results of particle-in-cell numerical simulations. We 
used the numerical code OSIRIS, described elsewhere 
l)Fonseca. et all l2002|h and we have performed 2D sim- 
ulations (1280 x 1280 cells, 128.0 x 128.0 (c/u> p , e ) 2 , 9 
particles/ (cell x species), 4 species) of the collision of 
weakly and fully relativistic neutral shells (electron- 
positron - e~e + , and electron-proton - e~p) moving in the 
z direction, across the xy simulation plan e, with parame- 
ters similar to those in ISilva. e t al. (2003). The temporal 
evolution of As as measured in the simulations is shown 
in Figure |3 Both a non-power-law nonrelativistic regime 
(until t « 10 - 20/w Pi e) and a power-law regime are clearly 
seen. The power-law fits yield Ab(£) oc t a with a s» 0.7 
for the sub-relativistic cases (it = 0.6c) and a w 1.1 for 
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the relativistic cases (7 = 10). Note also that the sec- 
ond power-law segment with the same index is present at 
t ;> 100/cjp,e in e~p sub-relativistic run, indicating pro- 
ton filament coalescences. A si milar behavior wa s also 
observed in 3D simulations (cf. ISilva. et "all 120 031 . but 
the significantly larger simulation planes presented here 
allow for improved statistics. At late times t ;> lOO/w^ei 
the evolution of A^ rolls off and is slower: a return cur- 
rent is set-up around each filament, decreasing the range 
of the magnetic field to just a few electron collisionless 
skin depths (the thickness of the region where the return 
current is flowing), thus partially shielding the filaments. 
Filament coalescence then occurs at a slower rate. 

5. DISCUSSION 

The two-stream instability generates magnetic fields at 
shock fronts very fast, with the typical e-folding time 
^growth ~ ~ 10 _5 (r/n) 1 / 2 s, where n is the particle 
(e.g., electron) number density in cm -3 , which is of order 
of 10~ 10 s for internal GRB shocks and is of order of 10 -4 s 
for external shocks. The characteristic spatial scale above 
which the field is essentially random, as predicted by the 
linear instability theory, is Aii n ~ cT growt h , which is of or- 
der of 10 cm and 10 7 cm for internal and external shocks, 
respectively. The extremely short spatial scales, i.e., sharp 
field gradients, must be rapidly destroyed by dissipation. 
Indeed, it would be the natural result of pitch-angle dif- 
fusion of current-carrying charges in the chaotic magnetic 
fields. So, the question arises: Why the generated mag- 
netic fields do not rapidly decay back to zero as soon as the 
instability shuts off? The answer is: The produced fields 
and the corresponding currents sel f-organize and for m a 
quasi- two-dimensional distribution (Sil va. et al.ll2003|) . A 
typical magnetic field gradient scale grows with time very 
rapidly, with approximately the light crossing time oc t; 
whereas the particle diffusion is a substantially slower pro- 
cess. Hence, diffusive dissipation is drastically reduced. To 
illustrate this, we consider the field diffusion equation 

d t B = ~nd 2 xx B (26) 



with the dissipation coefficient, k, being constant, for sim- 
plicity. Approximating the spatial derivative as d x ~ 
As(t) -1 ~ Aq (t/to)~ a , where Ao and to are constants 
and a > 1/2, we obtain the scaling of B with time as 

B(t) = B exp {{t 1 - 2 " - tl- 2a )/T l - 2a ) -> const., (27) 

as t -> 00, where t 1 " 2 " = (2a - 1)A§/(k^ q ). 

We also note that in some respect, the field scale 
growth is analogous to the inverse cascade in two- 
dimension al magnetohydrodynamic (MHD) turbulence 
(see, e.g., iBiskamp fc Bremerlll994f) . The crucial differ- 
ence is, however, the entirely kinetic nature of the process; 
at such small scales ~ c/lo p the MHD approximation is 
completely inapplicable. 

To conclude, in this paper we analytically investigated 
the long-time nonlinear dynamics of magnetic fields cre- 
ated at collisionless shocks by the two-stream instability. 
We demonstrated that the field correlation scale grows 
first exponentially and then nearly linearly with time, see 
Eqs. (|15f) and (|23|l . The transition from one regime 
to another occurs after few plasma times, see Eq. Q25[l. 
We compare our theoretical results with the state-of-the- 
art PIC numerica l simulation. Our fully 3D simulations 
ISilva. et al]l20f)! prove that the present simplified 2D 
analysis is accurate in the nonlinear regime until at least 
t ~ few x \Q 2 /oj p e . Whether (or when) our 2D model 
breaks down at later times cannot be tested with the 
present computer capabilities. The effect of the field evolu- 
tion on the observed spectrum and whether it can explain 
the spectral variability of the prompt GRB emission de- 
serves special consideration and will be discussed in the 
subsequent paper. 
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Fig. 1. — A schematic representation of the initial state. 

Fig. 2. — Temporal evolution of the field correlation length in four scenarios. The subscript R denotes relativistic shells (u^{ = 7V = ±10.0 cz 
while u = ±0.6cz otherwise, with u th = "fv th = 0.1c). In the range 20/u> p ^ < t < 100/w p>e , oc t a with a = 0.7, 0.6, 1, 1.2 for e~e+, e~p 
(e~e+)fl, (e~p)n, respectively. In the inset, we show the evolution of the B-field cquipartition parameter; note that at 20/up ie < t < 100/w Pl e 
eg ~ const.. 



